#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use uniquify;
use DBI;
use DBD::Pg;
use Data::Dump qw(dump);

(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    get_gene_positions_of_orthologs

	Description:

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------




sub get_synteny($$$\%\%\%\%)
{
	my ($dbh, $ortholog_type, $protocol_id, $genepos1, $genepos2, 
		$gene_ids1, $gene_ids2) = @_;
	print STDERR "\t\tGet gene positions for $ortholog_type...\n";

	print STDERR "\t\tGet synteny :\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT
		gl1.gene_position,
		gl2.gene_position,
        gl1.gene_id,
        gl2.gene_id
    FROM
		orthologs.ortholog_sets os1 NATURAL JOIN
		taxon.gene_loci gl1,
		orthologs.ortholog_sets os2 NATURAL JOIN
		taxon.gene_loci gl2
	WHERE
		os1.ortholog_type = '$ortholog_type' and
		os2.ortholog_type = '$ortholog_type' and
        os1.ortholog_id = os2.ortholog_id and
        os1.protocol_id = os2.protocol_id and
        os1.protocol_id = $protocol_id and
		gl1.chromosome ~ '^[0-9XYZW]+\$' and
		gl2.chromosome ~ '^[0-9XYZW]+\$' AND 
        os1.species = '$ortho_name1' AND 
        os2.species = '$ortho_name2'
    ORDER BY 
		gl1.gene_position,
		gl2.gene_position
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
	for (@$data)
	{
		push(@{$genepos1->{$ortholog_type}}, $_->[0]);
		push(@{$genepos2->{$ortholog_type}}, $_->[1]);
		push(@{$gene_ids1->{$ortholog_type}}, $_->[2]);
		push(@{$gene_ids2->{$ortholog_type}}, $_->[3]);
	}
}


#_________________________________________________________________________________________
#
#	Get chromosome Loci
#
#
sub get_chromosome_loci($\@\@\@$$)
{
	print STDERR ("\t\tRetrieving chromosome start ends...\n");
	my ($dbh, $chromosomes, $min_pos, $max_pos, $species, $sql_filter) = @_;
	$sql_filter = 'AND '.$sql_filter if length($sql_filter);
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT * from
    (
        SELECT
    		chromosome,
    		MIN(gene_position) - 0.5 AS gene_pos1,
    		MAX(gene_position) + 0.5 AS gene_pos2
    	FROM
    		gene_loci NATURAL JOIN ens_id
    	WHERE
    		chromosome ~ '^[0-9XYZW]+\$' AND 
            species = '$species'
            $sql_filter
    	GROUP BY
    		chromosome
     ) AS a
     ORDER BY gene_pos1;
PL/SQLCMD
	my $data =  $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
	print STDERR "\t\tSaving\t", scalar @$data, " chromosome loci...\n";
	for my $row(@$data)
	{
		push(@$chromosomes, $row->[0]);
		push(@$min_pos,  	$row->[1]);
		push(@$max_pos,  	$row->[2]);
	}
}



sub min(@)
{
	my @sorted = sort {$a<=>$b} @_;
	return $sorted[0];
}

sub max(@)
{
	my @sorted = sort {$a<=>$b} @_;
	return $sorted[-1];
}

sub filter_and_map_xy($$$$$$$$)
{
	my ($x, $y,  
		$x_beg, $x_end, 
		$y_beg, $y_end,
		$mapping_factor_x, $mapping_factor_y) = @_;

	die unless (scalar @$x == scalar @$y);

	my @new_x;
	my @new_y;
	for (my $i = 0; $i < scalar @$x; ++$i)
	{
		next if (	$x->[$i] < $x_beg || 
					$x->[$i] > $x_end ||
					$y->[$i] < $y_beg || 
					$y->[$i] > $y_end );
		
		push(@new_x, ($x->[$i] - $x_beg) * $mapping_factor_x);
		push(@new_y, ($y->[$i] - $y_beg) * $mapping_factor_y);
	}
	@$x = @new_x;
	@$y = @new_y;
}

sub filter_and_map(\@$$)
{
	my ($x, $beg, $mapping_factor) = @_;
	my @new_x = map {($_ - $beg) * $mapping_factor} @$x;
	@$x = @new_x;
}



#_________________________________________________________________________________________
#
#	print out chromosome grid
#
#_________________________________________________________________________________________
sub print_chromosome_grid(\@\@$)
{
	my ($x_chrm_loci, $y_chrm_loci, $directory) = @_;

	system ("mkdir -p $dir_orthologs_output/gene_synteny_map/$directory");
	open CHRMX, ">$dir_orthologs_output/gene_synteny_map/$directory/chromosomes_grid.txt" or die;
	
	for my $x_chrm_loci (@$x_chrm_loci)
	{
		print CHRMX "$x_chrm_loci\t0\t$x_chrm_loci\t1.0\r\n";
	}
	print CHRMX "0\t0\t0\t1.0\r\n";


	for my $y_chrm_loci (@$y_chrm_loci)
	{
		print CHRMX "1.0\t$y_chrm_loci\t0\t$y_chrm_loci\r\n";
	}
	print CHRMX "1.0\t0\t0\t0\r\n";
}


#_________________________________________________________________________________________
#
#	print out gene positions
#
#_________________________________________________________________________________________
sub print_gene_positions(\%\%\%\%$)
{
	my ($gene_pos_x, $gene_pos_y,
		$gene_id_x,$gene_id_y,
		$directory) = @_;
	
	system ("mkdir -p $dir_orthologs_output/gene_synteny_map/$directory");

	#
	#	print out genes by ortholog types
	#
	for my $ortholog_type('1 to 1', '1 to m', 'm to 1', 'm to m')
	{
		my $x = $gene_pos_x->{$ortholog_type};
		my $y = $gene_pos_y->{$ortholog_type};
		my $x_id = $gene_id_x->{$ortholog_type};
		my $y_id = $gene_id_y->{$ortholog_type};

		{
			open DOTS, ">$dir_orthologs_output/gene_synteny_map/$directory/$ortholog_type.txt" or die;
			for my $i(0.. (scalar @$x - 1))
			{
				print DOTS "$x->[$i]\t$y->[$i]\r\n";
			}
		}
		{
			open DOTS, ">$dir_orthologs_output/gene_synteny_map/$directory/$ortholog_type.gene_id.txt" or die;
			for my $i(0.. (scalar @$x - 1))
			{
				print DOTS	"$x->[$i]\t$y->[$i]\t".
							"$x_id->[$i]\t$y_id->[$i]\n";
			}
		}
	};
}


#_________________________________________________________________________________________
#
#	print out chromosome names
#
#_________________________________________________________________________________________
sub print_chromosome_names(\@\@\@\@\@\@$)
{

	my ($x_chrm_beg,  
		$x_chrm_end,  
		$x_chrm_names,
		$y_chrm_beg,  
		$y_chrm_end,
		$y_chrm_names,
		$directory) = @_;

	system ("mkdir -p $dir_orthologs_output/gene_synteny_map/$directory");

	open CHRMX, ">$dir_orthologs_output/gene_synteny_map/$directory/x_chromosome_names.txt" or die;
	for my $i (0.. $#$x_chrm_names)
	{
		my $pos = ($x_chrm_beg->[$i] + $x_chrm_end->[$i]) / 2; 
		print CHRMX "$pos\t-0.05\t$x_chrm_names->[$i]\r\n";
	}

	open CHRMY, ">$dir_orthologs_output/gene_synteny_map/$directory/y_chromosome_names.txt" or die;
	for my $i (0.. $#$y_chrm_names)
	{
		my $pos = ($y_chrm_beg->[$i] + $y_chrm_end->[$i]) / 2; 
		print CHRMY "-0.05\t$pos\t$y_chrm_names->[$i]\r\n";
	}
}


sub map_array(\@\%)
{
	my ($array, $map) =@_;
	my @new_array = map {$map->{$_}} @$array;
	@$array = @new_array;
}


#_________________________________________________________________________________________
#
#	collage_and_print_syntenic_map
#
#_________________________________________________________________________________________
sub collage_and_print_syntenic_map($$$$)
{
	my ($dbh, $protocol_id, $directory, $sql_filter) = @_;
	print STDERR "\tSynteny for $directory...\n";

	#
	#	Get Loci (chromosomes and genes)
	#
	
		#
		#	get chromosomes	start and end
		#
		my @x_chrm_names;
		my @y_chrm_names;
		my (@x_chrm_beg, @x_chrm_end);
		my (@y_chrm_beg, @y_chrm_end);
		get_chromosome_loci($dbh, @x_chrm_names, @x_chrm_beg, @x_chrm_end, $ortho_name1, $sql_filter);
		get_chromosome_loci($dbh, @y_chrm_names, @y_chrm_beg, @y_chrm_end, $ortho_name2, $sql_filter);
	
		# make sure start and end chromosome lines are not drawn twice
		my @x_chrm_loci = sort {$a <=> $b} uniquify (@x_chrm_beg, @x_chrm_end);
		my @y_chrm_loci = sort {$a <=> $b} uniquify (@y_chrm_beg, @y_chrm_end);
		
		
		#
		#	get synteny for each ortholog type
		#
		my %gene_pos_x;
		my %gene_pos_y;
		my %gene_id_x;
		my %gene_id_y;
		for my $ortholog_type('1 to 1', '1 to m', 'm to 1', 'm to m')
		{
			get_synteny($dbh,  $ortholog_type, $protocol_id, %gene_pos_x, %gene_pos_y, 
						%gene_id_x, %gene_id_y);
		}


	#
	#	map gene positions to successive indices (sans gaps)
	#		i.e. 1.5, 1.6, 2.8 => 1, 2, 3 etc.
	#
	my (%old_x_to_new, %old_y_to_new);
	{	
		#
		#	put all gene positions and chromosomes into one big array
		#
		my @all_loci_x;
		my @all_loci_y;
		for my $ortholog_type('1 to 1', '1 to m', 'm to 1', 'm to m')
		{
			push(@all_loci_x, @{$gene_pos_x{$ortholog_type}});
			push(@all_loci_y, @{$gene_pos_y{$ortholog_type}});
		}
		die unless scalar @all_loci_x == @all_loci_y;
		push(@all_loci_x, @x_chrm_loci);
		push(@all_loci_y, @y_chrm_loci);
	
	
		
		#
		#	use one big array to map gene positions to successive indices sans gaps
		#
		@all_loci_x = sort{$a <=> $b} uniquify(@all_loci_x);
		@all_loci_y = sort{$a <=> $b} uniquify(@all_loci_y);
		@old_x_to_new{@all_loci_x} = (0..$#all_loci_x);
		@old_y_to_new{@all_loci_y} = (0..$#all_loci_y);
	}	




	#
	#	Convert gene positions to gapless positions
	#
	{
		# chromosome loci
		map_array(@x_chrm_beg , %old_x_to_new);
		map_array(@x_chrm_end , %old_x_to_new);
		map_array(@x_chrm_loci, %old_x_to_new);
		map_array(@y_chrm_beg , %old_y_to_new);
		map_array(@y_chrm_end , %old_y_to_new);
		map_array(@y_chrm_loci, %old_y_to_new);

		# map for orthologues
		for my $ortholog_type('1 to 1', '1 to m', 'm to 1', 'm to m')
		{
			map_array(@{$gene_pos_x{$ortholog_type}} , %old_x_to_new);
			map_array(@{$gene_pos_y{$ortholog_type}} , %old_y_to_new);
		}
	}
	#print "x_chrm_beg \n", dump (@x_chrm_beg ), "\n\n";
	#print "x_chrm_end \n", dump (@x_chrm_end ), "\n\n";
	#print "x_chrm_loci\n", dump (@x_chrm_loci), "\n\n";
	#print "y_chrm_beg \n", dump (@y_chrm_beg ), "\n\n";
	#print "y_chrm_end \n", dump (@y_chrm_end ), "\n\n";
	#print "y_chrm_loci\n", dump (@y_chrm_loci), "\n\n";
	#print "gene_pos_x \n", dump (%gene_pos_x ), "\n\n";
	#print "gene_pos_y \n", dump (%gene_pos_y ), "\n\n";

	#
	#	prepare mapping to 0.0 <-> 1.0
	#

		# get minimum and maximum position of axes
		my ($x_beg) = min(@x_chrm_beg);
		my ($x_end) = max(@x_chrm_end);
		my ($y_beg) = min(@y_chrm_beg);
		my ($y_end) = max(@y_chrm_end);
	
		# now convert all to 0..1 using positions from axis min end
		my $mapping_factor_x = 1.0 / ($x_end - $x_beg);
		my $mapping_factor_y = 1.0 / ($y_end - $y_beg);


	#
	#	map to 0.0 <-> 1.0
	#

		# only keep gene positions which will fall between 0 and 1
		for my $ortholog_type('1 to 1', '1 to m', 'm to 1', 'm to m')
		{
			filter_and_map_xy(	$gene_pos_x{$ortholog_type}, 
								$gene_pos_y{$ortholog_type},
								$x_beg,
								$x_end,
								$y_beg,
								$y_end,
								$mapping_factor_x,
								$mapping_factor_y);
		};

		# convert all chromosome positions
		filter_and_map(@x_chrm_beg, $x_beg, $mapping_factor_x);
		filter_and_map(@x_chrm_end, $x_beg, $mapping_factor_x);
		filter_and_map(@x_chrm_loci, $x_beg, $mapping_factor_x);
	
		filter_and_map(@y_chrm_beg, $y_beg, $mapping_factor_y);
		filter_and_map(@y_chrm_end, $y_beg, $mapping_factor_y);
		filter_and_map(@y_chrm_loci, $y_beg, $mapping_factor_y);
	

		#print ("*" x 80), "\n\nFILTER 0->1.0\n\n";
		#print "x_chrm_beg \n", dump (@x_chrm_beg ), "\n\n";
		#print "x_chrm_end \n", dump (@x_chrm_end ), "\n\n";
		#print "x_chrm_loci\n", dump (@x_chrm_loci), "\n\n";
		#print "y_chrm_beg \n", dump (@y_chrm_beg ), "\n\n";
		#print "y_chrm_end \n", dump (@y_chrm_end ), "\n\n";
		#print "y_chrm_loci\n", dump (@y_chrm_loci), "\n\n";
		#print "gene_pos_x \n", dump (%gene_pos_x ), "\n\n";
		#print "gene_pos_y \n", dump (%gene_pos_y ), "\n\n";



	#
	#	print
	#
		print_gene_positions(%gene_pos_x, %gene_pos_y, %gene_id_x, %gene_id_y, $directory);
		print_chromosome_grid(@x_chrm_loci, @y_chrm_loci, $directory);
		print_chromosome_names(@x_chrm_beg, @x_chrm_end, @x_chrm_names, 
							   @y_chrm_beg, @y_chrm_end, @y_chrm_names,
							   $directory);
	
}


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Get data from panda
#

my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;


collage_and_print_syntenic_map($dbh, $curr_protocol_id, "all", "");
#print_synteny($dbh, $curr_protocol_id, "M4HX",	"(species = '$ortho_name1' and chromosome = 'X') or".
#												"(species = '$ortho_name2' and chromosome = '4')");
#print_synteny($dbh, $curr_protocol_id, "M7HX",	"(species = '$ortho_name1' and chromosome = 'X') or".
#												"(species = '$ortho_name2' and chromosome = '7')");

print STDERR "Finished.\n";

